# Project: Serbia Survey
# Authors: William O'Brochta and Patrick Cunha Silva
# Goal: Analysis on voters' perceptions of politicians (H2)

# Load Data
load(file = here('VotersTowardPoliticians.RData'))

##########################################################
################# Main Models for H2 #####################
##########################################################

#------------------ Unified DV ------------------------#
# Model for H2, no controls
m1a <- lm(Y ~ slogan_cyrillic, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))
ub1a <- qt(0.975, df = df.residual(m1a)) * se1a['slogan_cyrillic'] + coef(m1a)['slogan_cyrillic']
lb1a <- qt(0.025, df = df.residual(m1a)) * se1a['slogan_cyrillic'] + coef(m1a)['slogan_cyrillic']

#------------------ better_world ------------------------#
# Model for H2, no controls
m1b <- lm(better_world ~ slogan_cyrillic, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))
ub1b <- qt(0.975, df = df.residual(m1b)) * se1b['slogan_cyrillic'] + coef(m1b)['slogan_cyrillic']
lb1b <- qt(0.025, df = df.residual(m1b)) * se1b['slogan_cyrillic'] + coef(m1b)['slogan_cyrillic']

#------------------ most_places ------------------------#
# Model for H2, no controls
m1c <- lm(most_places ~ slogan_cyrillic, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))
ub1c <- qt(0.975, df = df.residual(m1c)) * se1c['slogan_cyrillic'] + coef(m1c)['slogan_cyrillic']
lb1c <- qt(0.025, df = df.residual(m1c)) * se1c['slogan_cyrillic'] + coef(m1c)['slogan_cyrillic']

#------------------ influence ------------------------#
# Model for H2, no controls
m1d <- lm(influence ~ slogan_cyrillic, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))
ub1d <- qt(0.975, df = df.residual(m1d)) * se1d['slogan_cyrillic'] + coef(m1d)['slogan_cyrillic']
lb1d <- qt(0.025, df = df.residual(m1d)) * se1d['slogan_cyrillic'] + coef(m1d)['slogan_cyrillic']


# Effect size
coef(m1b)['slogan_cyrillic']/mean(mydf$better_world)
coef(m1c)['slogan_cyrillic']/mean(mydf$most_places)
coef(m1d)['slogan_cyrillic']/mean(mydf$influence)
coef(m1a)['slogan_cyrillic']/mean(mydf$Y)

coef(m1b)['slogan_cyrillic']/sd(mydf$better_world)
coef(m1c)['slogan_cyrillic']/sd(mydf$most_places)
coef(m1d)['slogan_cyrillic']/sd(mydf$influence)
coef(m1a)['slogan_cyrillic']/sd(mydf$Y)

# Plot
par(mar = c(5, 12, 2, 2))
plot(y = c(0.75, 2, 3.25, 4.5), x = c(coef(m1b)['slogan_cyrillic'], 
                                      coef(m1c)['slogan_cyrillic'], coef(m1d)['slogan_cyrillic'], coef(m1a)['slogan_cyrillic']), 
     xlim = c(-0.6, 0.6), ylim = c(0, 5),
     pch = 20, cex = 2.25, 
     ylab = '', xlab = 'ATE of Cyrillic Script', cex.lab = 1.75, 
     axes = FALSE)
lines(y = c(0.75, 0.75), x = c(ub1b, lb1b), lwd = 1.5)
lines(y = c(2, 2), x = c(ub1c, lb1c), lwd = 1.5)
lines(y = c(3.25, 3.25), x = c(ub1d, lb1d), lwd = 1.5)
lines(y = c(4.5, 4.5), x = c(ub1a, lb1a), lwd = 1.5)
axis(1, las = 1, cex.axis = 1.25)
axis(2, labels = c('Better World\nif Equal to Serbia', ' Serbia is Better\nThan Most Places', 'More Serbian\nInfluence is Better', 'Nationalism\nIndex'), 
     at = c(0.75, 2, 3.25, 4.5), las = 2, tick = FALSE, cex.axis = 1.5)
abline(v = 0, lty = 2)


###########################################################################
################# SI A.1 - Descriptive Stats - Part 2 #####################
###########################################################################

# Create dummies
mydf <- dummy_cols(mydf, select_columns = c('gender', 'ageCat', 'other_languages'))


# List of covariates
covariates <- c("cyrillic_second",
                "cyrillic_start",
                'milena',
                "gender_woman",
                "gender_man",
                "gender_no_response",
                "college",
                "ageCat_18-35",
                "ageCat_36-49",
                "ageCat_50-64",
                "ageCat_65+",
                "married",
                "unemployed",
                "urban",
                "other_languages_0",
                "other_languages_1",
                "other_languages_2",
                "other_languages_>2",
                "dif_lc_home",
                "dif_lc_work",
                'cyrillicImportance',
                "serbian",
                "sns",
                "polKnowQ1",
                "polKnowQ2",
                "intPolitics")

# Final List of variables
variables <- c("Y", "better_world", "most_places", "influence", "slogan_cyrillic", covariates)

# List of labels
var.names <- c('Nationalism Index (Politician)',
               'More Serbian influence is Better (Politician)', 
               'Serbia is Better than Most Places (Politician)',
               'World would be Better if it was Equal to Serbia (Politician)',
               'Second Randomization (Cyrillic Slogan)',
               'First Randomization (Cyrillic)',
               'Alphabet Ad (Cyrillic)', 
               'Politician = Milena',
               "Gender = Woman",
               "Gender = Man", 
               "Gender = Prefer not to share", 
               "Education = College", 
               "Age = 18-35", 
               "Age = 36-49", 
               "Age = 50-64", 
               "Age = 65+", 
               "Marital status = Married", 
               "Unemployed", 
               "Urban", 
               "# of languages spoken = 0",
               "# of languages spoken = 1", 
               "# of languages spoken = 2", 
               "# of languages spoken > 2", 
               "Use of Cyrillic-Latin (at home)", 
               "Use of Cyrillic-Latin (at work)", 
               "Importance of Cyrillic",
               "Serbian", 
               "SNS",
               "Knows the # of MPs", 
               "Knows the presidential term length", 
               "Interest in politics")

# Generate Table
stargazer::stargazer(mydf[, variables], covariate.labels = var.names, type = 'text')

#########################################################################
################# SI A.11 - H2: Randomization Checks ####################
#########################################################################

# List of covariates
covariates <- c("cyrillic_second",
                "cyrillic_start",
                'milena',
                "gender_woman",
                "gender_man",
                "gender_no_response",
                "college",
                "ageCat_18-35",
                "ageCat_36-49",
                "ageCat_50-64",
                "ageCat_65+",
                "married",
                "unemployed",
                "urban",
                "other_languages_0",
                "other_languages_1",
                "other_languages_2",
                "other_languages_>2",
                "dif_lc_home",
                "dif_lc_work",
                'cyrillicImportance',
                "serbian",
                "sns",
                "polKnowQ1",
                "polKnowQ2",
                "intPolitics")

# List of variable names
var.names <- c('First Randomization (Cyrillic)',
               'Alphabet Ad (Cyrillic)', 
               'Politician = Milena',
               "Gender = Woman",
               "Gender = Man", 
               "Gender = Prefer not to share", 
               "Education = College", 
               "Age = 18-35", 
               "Age = 36-49", 
               "Age = 50-64", 
               "Age = 65+", 
               "Marital status = Married", 
               "Unemployed", 
               "Urban", 
               "# of languages spoken = 0",
               "# of languages spoken = 1", 
               "# of languages spoken = 2", 
               "# of languages spoken > 2", 
               "Use of Cyrillic-Latin (at home)", 
               "Use of Cyrillic-Latin (at work)", 
               "Importance of Cyrillic",
               "Serbian", 
               "SNS",
               "Knows the # of MPs", 
               "Knows the presidential term length", 
               "Interest in politics")

# Save pvalues
pvalues <- list()
for(i in 1:length(covariates)){pvalues[[i]] <- t.test(mydf[, covariates[i]] ~ mydf$slogan_cyrillic)$p.value}
unlist(pvalues)

par(mar = c(6, 15, 1, 1))
plot(x = pvalues, y = 1:length(covariates), 
     xlim = c(0, 1), 
     ylim = c(1,length(covariates)),
     pch = 20, 
     xlab = 'p-value',
     ylab = '', 
     axes = FALSE, 
     cex = 1.75, 
     cex.lab = 1.5)
axis(1, at = seq(0, 1, by = 0.05))
axis(2, at = 1:length(covariates), labels = var.names, las = 1)
abline(v = 0.01, lty = 3)
abline(v = 0.05, lty = 3)


###########################################################################
################# SI A.12 - H2: Complete Results ##########################
###########################################################################

# Table for APSR
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     style = 'apsr', omit = 'Constant', type = 'text')


#######################################################################
############ SI A.13 - H2: Descriptive Figures (Outcomes) ##############
#######################################################################


######### "The world would be a better place if people from other countries were more like Serbians” 
layout(matrix(c(1), nrow = 1, ncol = 1))
t0 <- unname(table(mydf$better_world[mydf$slogan_cyrillic == 0]))
t1 <- unname(table(mydf$better_world[mydf$slogan_cyrillic == 1]))
mylabels <- c('Strongly Disagree', 'Disagree', "Doesn't Agree or Disagree", 'Agree', 'Strongly Agree')
par(mar = c(15, 5, 2, 2))
m1 <- barplot(prop.table(t1) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m1, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m0 <- barplot(prop.table(t0) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m10 <- barplot((prop.table(t1) - prop.table(t0)) * 100, ylim = c(-10, 10), ylab = 'Difference in % of Respondents',  border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -10, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)

######### "Generally speaking, Serbia is no better than most other countries"
t0 <- unname(table(mydf$most_places[mydf$slogan_cyrillic == 0]))
t1 <- unname(table(mydf$most_places[mydf$slogan_cyrillic == 1]))
mylabels <- c('Strongly Disagree', 'Disagree', "Doesn't Agree or Disagree", 'Agree', 'Strongly Agree')
par(mar = c(15, 5, 2, 2))
m1 <- barplot(prop.table(t1) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m1, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m0 <- barplot(prop.table(t0) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m10 <- barplot((prop.table(t1) - prop.table(t0)) * 100, ylim = c(-10, 10), ylab = 'Difference in % of Respondents',  border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -10, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)

######### "Generally, the more influence Serbia has on other nations, the better off they are"
t0 <- unname(table(mydf$influence[mydf$slogan_cyrillic == 0]))
t1 <- unname(table(mydf$influence[mydf$slogan_cyrillic == 1]))
mylabels <- c('Strongly Disagree', 'Disagree', "Doesn't Agree or Disagree", 'Agree', 'Strongly Agree')
par(mar = c(15, 5, 2, 2))
m1 <- barplot(prop.table(t1) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m1, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m0 <- barplot(prop.table(t0) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m10 <- barplot((prop.table(t1) - prop.table(t0)) * 100, ylim = c(-10, 10), ylab = 'Difference in % of Respondents',  border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -10, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)

###########################################################################
################# SI A.14 - H2: Models with Covariates #####################
###########################################################################

#------------------ Unified DV ------------------------#
# Model for H2, with controls
m1a <- lm(Y ~ slogan_cyrillic + cyrillic_second + cyrillic_start  + milena +
            gender + college + ageCat + married + unemployed + urban +
            other_languages + dif_lc_home + dif_lc_work + cyrillicImportance +
            serbian + sns + polKnowQ1 + polKnowQ2 + intPolitics, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H2, with controls
m1b <- lm(better_world ~ slogan_cyrillic + cyrillic_second + cyrillic_start  + milena +
            gender + college + ageCat + married + unemployed + urban +
            other_languages + dif_lc_home + dif_lc_work + cyrillicImportance +
            serbian + sns + polKnowQ1 + polKnowQ2 + intPolitics, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H2, with controls
m1c <- lm(most_places ~ slogan_cyrillic + cyrillic_second + cyrillic_start  + milena +
            gender + college + ageCat + married + unemployed + urban +
            other_languages + dif_lc_home + dif_lc_work + cyrillicImportance +
            serbian + sns + polKnowQ1 + polKnowQ2 + intPolitics, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H2, with controls
m1d <- lm(influence ~ slogan_cyrillic + cyrillic_second + cyrillic_start  + milena +
            gender + college + ageCat + married + unemployed + urban +
            other_languages + dif_lc_home + dif_lc_work + cyrillicImportance +
            serbian + sns + polKnowQ1 + polKnowQ2 + intPolitics, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

# Generate Table
stargazer(m1a, m1b, m1c, m1d, 
          se = list(se1a, se1b, se1c, se1d), 
          type = 'text',
          covariate.labels = c('Treatment 2 (Slogan in Cyrillic)', 
                               'Treatment 1 (Survey in Cyrillic)',
                               'Alphabet Ad (Cyrillic)', 
                               "Politician = Milena",
                               "Gender = Man", 
                               "Gender = Prefer not to share", 
                               "Education = College", 
                               "Age = 36-49", 
                               "Age = 50-64", 
                               "Age = 65+", 
                               "Marital Status = Married", 
                               "Unemployed", 
                               "Urban", 
                               "# of languages spoken = 1", 
                               "# of languages spoken = 2", 
                               "# of languages spoken > 2", 
                               "Delta use of Cyrillic to Latin (at home)", 
                               "Delta use of Cyrillic to Latin (at work)", 
                               "Importance of Cyrillic",
                               "Ethnicity = Serbian", 
                               "Party Identification = SNS",
                               "Knows the # of MPs", 
                               "Knows how long is the presidential term", 
                               "Interest in Politics"),
          star.cutoffs = c(0.05, 0.01),
          omit.table.layout = "n",	
          dep.var.labels = c('Nationalism index', "More Serbian is better", 
                             "Serbia is bettern than most places", "Better World if Equal to Serbia"),
          title = "Effect of Cyrillic Alphabet on Respondents' Perceptions of Politicians, Including Covariates", 
          single.row = TRUE,
          no.space = TRUE ,
          omit.stat = c("ser", 'f'),
          style = 'apsr'
)



############################################################################
################# SI A.15 - H2: Result by Ad's Alphabet ####################
############################################################################

# Create Dummies (start_treat_slogan)
mydf$latin_latin <- as.numeric(mydf$latin_start == TRUE & mydf$slogan_latin == TRUE)
mydf$latin_cyrillic <- as.numeric(mydf$latin_start == TRUE & mydf$slogan_latin == FALSE)
mydf$cyrillic_latin <- as.numeric(mydf$latin_start == FALSE & mydf$slogan_latin == TRUE)
mydf$cyrillic_cyrillic <- as.numeric(mydf$latin_start == FALSE & mydf$slogan_latin == FALSE)

# Latin + Latin is the control

#------------------ Unified DV ------------------------#
# Model for H2, no controls
m1a <- lm(Y ~ latin_cyrillic + cyrillic_latin + cyrillic_cyrillic, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H2, no controls
m1b <- lm(better_world ~ latin_cyrillic + cyrillic_latin + cyrillic_cyrillic, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H2, no controls
m1c <- lm(most_places ~ latin_cyrillic + cyrillic_latin + cyrillic_cyrillic, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H2, no controls
m1d <- lm(influence ~ latin_cyrillic + cyrillic_latin + cyrillic_cyrillic, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

# Table for APSR
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     covariate.labels = c('Ad = Latin, Slogan = Cyrillic',
                                          'Ad = Cyrillic, Slogan = Latin',
                                          'Ad = Cyrillic, Slogan = Cyrillic'),
                     title = "Effect of Alphabet on Respondents' Perceptions of Politicians",
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr',
                     star.cutoffs = c(0.05, 0.01),
                     type = 'text'
)


############################################################################################################
############ SI A.16 - H2: Interaction Political Knowledge, Interest in Politics, and Party ID #############
############################################################################################################

### Interest in Politics

#------------------ Unified DV ------------------------#
# Model for H1, interaction with interest in politics
m1a <- lm(Y ~ slogan_cyrillic * intPolitics, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H1, interaction with interest in politics
m1b <- lm(better_world ~ slogan_cyrillic * intPolitics, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H1, interaction with interest in politics
m1c <- lm(most_places ~ slogan_cyrillic * intPolitics, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H1, interaction with interest in politics
m1d <- lm(influence ~ slogan_cyrillic * intPolitics, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

# Table 
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     type = 'text',
                     covariate.labels = c('Slogan = Cyrillic',
                                          'Interest in Politics',
                                          'Slogan = Cyrillic x Interest in Politics'),
                     dep.var.labels = c('Nationalism index', "More Serbian is better", 
                                        "Serbia is better than most places", "Better World if Equal to Serbia"),                     
                     title = "Effect of Alphabet on Respondents' Perceptions of Politicians Conditional upon Interest in Politics",
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr',
                     star.cutoffs = c(0.05, 0.01)
)
## Create plots
model_list <- list(m1a, m1b, m1c, m1d)
sapply(model_list, plot_interaction, x = 'slogan_cyrillic', z = 'slogan_cyrillic:intPolitics', x_name = 'Cyrillic Slogan')

#### Political Knowledge

#------------------ Unified DV ------------------------#
# Model for H1, interaction with political knowledge
m1a <- lm(Y ~ slogan_cyrillic * polKnow, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H1, interaction with political knowledge
m1b <- lm(better_world ~ slogan_cyrillic * polKnow, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H1, interaction with political knowledge
m1c <- lm(most_places ~ slogan_cyrillic * polKnow, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H1, interaction with political knowledge
m1d <- lm(influence ~ slogan_cyrillic * polKnow, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

# Table 
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     type = 'text',
                     covariate.labels = c('Slogan = Cyrillic',
                                          'Political Knowledge',
                                          'Slogan = Cyrillic x Political Knowledge'),
                     dep.var.labels = c('Nationalism index', "More Serbian is better", 
                                        "Serbia is better than most places", "Better world if equal to Serbia"),                     
                     title = "Effect of Alphabet on Respondents' Perceptions of Politicians Conditional upon Political Knowledge",
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr',
                     star.cutoffs = c(0.05, 0.01)
)

## Create plots
model_list <- list(m1a, m1b, m1c, m1d)
sapply(model_list, plot_interaction, z_values = 0:2, x = 'slogan_cyrillic', z = 'slogan_cyrillic:polKnow', z_name = 'Political Knowledge', x_name = 'Cyrillic Slogan')

#### Party ID

#------------------ Unified DV ------------------------#
# Model for H1, interaction with partisan preference
m1a <- lm(Y ~ slogan_cyrillic * sns, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H1, interaction with  partisan preference
m1b <- lm(better_world ~ slogan_cyrillic * sns, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H1, interaction with partisan preference
m1c <- lm(most_places ~ slogan_cyrillic * sns, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H1, interaction with partisan preference
m1d <- lm(influence ~ slogan_cyrillic * sns, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

# Table 
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     type = 'text',
                     covariate.labels = c('Slogan = Cyrillic',
                                          'Party Identification = SNS',
                                          'Slogan = Cyrillic x Party Identification = SNS'),
                     dep.var.labels = c('Nationalism index', "More Serbian is better", 
                                        "Serbia is better than most places", "Better world if equal to Serbia"),                     
                     title = "Effect of Alphabet on Respondents' Perceptions of Politicians Conditional upon Party Identification = SNS",
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr',
                     star.cutoffs = c(0.05, 0.01)
)

## Create plots
model_list <- list(m1a, m1b, m1c, m1d)
sapply(model_list, plot_interaction, z_values = 0:1, x = 'slogan_cyrillic', z = 'slogan_cyrillic:sns', z_name = 'Party Identification = SNS', x_name = 'Cyrillic Slogan')
